/*
  Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes related to emission of rays.

#ifndef __emission__
#define __emission__

#include <algorithm> 
#include <vector> 
#include "mcrx-types.h"
#include "ray.h"
#include "sphparticle.h"
#include "montecarlo.h"
#include "random.h"
#include "constants.h"
#include "boost/shared_ptr.hpp"
#include "vecops.h"
#include "angular.h"
#include "grid.h"
#include <memory>
#include "almost_equal.h"
#include "chromatic_policy.h"
#include "optical.h"

namespace mcrx {
  template<typename, typename> class emission;
  template<typename, typename> class isotropic_emission;
  template<typename, template<typename> class, typename> 
  class emission_collection;
  template <typename, typename> class volumetric_emission;  
  template <typename T> class volumetric_emission<polychromatic_policy, T>;  

  template <typename, typename> class grid_cell_emission;
  template<typename, typename> class pointsource_emission;
  template<typename, typename> class biased_pointsource_emission;
  template<typename, typename> class laser_emission;
  template<typename, typename> class floodlight_emission;
  template<typename, typename> class divergent_floodlight_emission;
  template<typename, typename> class external_illumination;
  template<typename, typename> class particle_emission;
  template<typename, typename> class particles_emission;

  /** Traits class that defines the class for an emission distribution
      for different grids. */
  template<template <typename> class, typename, typename> class emitter_types {};

  /** Specialization of the emitter_types class for adaptive_grid. */
  template<typename> class adaptive_grid;
  template<typename chromatic_policy, 
	   typename rng_policy> class emitter_types<mcrx::adaptive_grid, 
						    chromatic_policy, 
						    rng_policy> {
  public:   
    typedef grid_cell_emission<chromatic_policy, rng_policy > T_emitter;
  };

#ifdef WITH_AREPO
  /** Specialization of the emitter_types class for arepo_grid. */
  template<typename> class arepo_grid;
  template<typename, typename> class arepo_cell_emission;
  template<typename chromatic_policy, 
	   typename rng_policy> class emitter_types<mcrx::arepo_grid, 
						    chromatic_policy, 
						    rng_policy> {
  public:   
    typedef arepo_cell_emission<chromatic_policy, rng_policy > T_emitter;
  };
#endif

  template <typename A, typename B>
  void assign(volumetric_emission<A,B>& lhs,
	      const volumetric_emission<A,B>& rhs);
  template <typename A>
  void assign(volumetric_emission<polychromatic_policy,A>& lhs,
	      const volumetric_emission<polychromatic_policy,A>& rhs);
  template <typename A, typename B>
  void assign(grid_cell_emission<A,B>& lhs,
	      const grid_cell_emission<A, B>& rhs);

}


/** Abstract base class for emission objects, which emit rays
    according to some design probability distribution.  Defines the
    interface.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::emission {
protected:

public:
  typedef typename chromatic_policy::T_lambda T_lambda;
  typedef typename chromatic_policy::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;

  virtual ~emission () {};
  /// Exception class, thrown if there is no emission in the grid.
  class no_emission {}; 
  
  /** The emission function, returns a ray and the angular
      distribution of the emitted rays.  The function takes the policy for
      generating random numbers as an argument, because the random number
      generator is thread-local.  */
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy&, T_float&) const=0; 

  virtual T_lambda zero_lambda() const=0;
  virtual T_float luminosity() const=0;
  virtual boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const=0;
};


/** Abstract base class for all emission classes with isotropic
    angular distribution.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::isotropic_emission : 
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;

protected:
  isotropic_angular_distribution<T_lambda> ad_;

public:
  isotropic_emission() {};
  isotropic_emission(const T_lambda& l): ad_(l) {};
};


/** The emission_collection class manages emission from a collection
    of emitters. The emitters themselves are not contained in this
    class, it only keeps pointers to them. For emission, it picks the
    individual emitters according to their statistical weight and then
    asks that emitter to emit the ray. */
template <typename chromatic_policy,  
	  template<typename> class sampling_policy, 
	  typename rng_policy_type>
class mcrx::emission_collection :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename mcrx::emission<chromatic_policy, rng_policy_type> T_emission;
  typedef typename T_emission::T_lambda T_lambda;
  typedef typename T_emission::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;  
  typedef ray<T_lambda> T_ray;
  typedef probability_distribution<T_emission, 
				   sampling_policy, 
				   T_rng_policy> T_distribution;

protected:
	    
  T_distribution distr_;
  /// Default constructor is only for derived classes.
  emission_collection () {};
  
public:
  emission_collection (const std::vector<emission<chromatic_policy, 
		       T_rng_policy>*>& emitters,
		       const array_1& weights) :
    emission<chromatic_policy, rng_policy_type>()
  { distr_.setup(emitters, weights); };

  /** The emit function simply defers to the individual
      emitters after drawing one. */
  virtual ray_distribution_pair<T_lambda> emit (T_biaser b,
						T_rng_policy& p,
						T_float& prob) const {
    assert(!distr_.pointers().empty());
    T_float temp_prob1, temp_prob2;
    ray_distribution_pair<T_lambda> 
      retval(distr_.sample(p, temp_prob1)->emit(b, p, temp_prob2));
    prob=temp_prob1*temp_prob2; return retval;
  };
  T_float total_emission_weight () const {
    return distr_.total_weight();};

  virtual T_float luminosity() const { return total_emission_weight(); };

  /** Returns a reference to the distribution object. */
  const T_distribution& distribution() const {return distr_;};
  
  /** This function asks the first emitter for its zero_lambda, it's
      dumb but easy. */
  T_lambda zero_lambda() const {
    assert(!distr_.pointers().empty());
    return distr_.pointers().front()->zero_lambda();};
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new emission_collection<chromatic_policy, sampling_policy, rng_policy_type>(*this));};
};


/** Emission class that handles isotropic emission from a finite sized
    particle. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::particle_emission : 
  public mcrx::isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef particle<T_lambda> T_particle;
  typedef rng_policy_type T_rng_policy;

private:
  /** The particle keeps the data about what is being emitted, like
      the position, velocity and emission of the emitter. */
  T_particle p; 

public:
  particle_emission (const vec3d& position, const vec3d& velocity, 
		     T_float size, T_lambda em):
    isotropic_emission<chromatic_policy, rng_policy_type>(em), 
    p (0, position, velocity, size, em) {};
  
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float&) const; 

  virtual T_float luminosity() const { return 1; };

  T_lambda zero_lambda() const {return T_lambda(0.*p.data());}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new particle_emission<chromatic_policy, rng_policy_type>(*this));};
};



/** Emission class representing emission from an isotropic point source. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::pointsource_emission : 
  public mcrx::isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  /// Position of the point source.
  vec3d position;

  T_lambda emission_; ///< The "thing" being emitted.
public:
  pointsource_emission (const vec3d& pos, const T_lambda em):
    isotropic_emission<chromatic_policy, rng_policy_type>(em),
    position (pos), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 

  virtual T_float luminosity() const { return 1; };

  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new pointsource_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** An isotropic point source which is biased towards the xy
    plane. This does not inherit from isotropic_emission because the
    probability distribution is different, even though it IS
    isotropic. Unfortunately we don't have a nice way of handling
    sources with biased emission, so this is a one-off case. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::biased_pointsource_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

private:
  vec3d position; ///< The position of the source.
  T_lambda emission_; ///< The "thing" being emitted.

public:
  biased_pointsource_emission (const vec3d& pos, 
			       const T_lambda  em ):
    position (pos), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new biased_pointsource_emission<chromatic_policy, rng_policy_type>(*this));};
};  


/** Emission class representing collimated emission from a point.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::laser_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  /// The angular distribution object.
  collimated_angular_distribution<T_lambda> d;
  
  T_lambda emission_; ///< The "thing" being emitted.

public:
  laser_emission (const vec3d& pos, const vec3d& dir, 
		  const T_lambda  em ):
    position (pos), d (dir, em), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new laser_emission<chromatic_policy, rng_policy_type>(*this));};
};




/** Emission class representing a finite-sized beam of collimated
    emission from a point, covering a certain areal cross-section
    (like an old aircraft searchlight).  */
template <typename chromatic_policy, typename rng_policy_type>
  class mcrx::floodlight_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
  public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  T_float radius; ///< The radius of the beam.
  collimated_angular_distribution<T_lambda> d;
  
  T_lambda emission_; // The "thing" being emitted.

public:
  floodlight_emission (const vec3d& pos, const vec3d& dir, T_float r, 
		       const T_lambda  em ):
    position (pos), d (dir, em), radius (r), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };
  
    virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new floodlight_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** Emission class representing an almost-collimated beam with
    specified (initial) cross section and divergence. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::divergent_floodlight_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  cone_angular_distribution<T_lambda> d_;
  T_float radius_; ///< The radius of the beam.
  T_float theta_; ///< The opening angle of the radiation.
  T_lambda emission_; // The "thing" being emitted.

public:
  divergent_floodlight_emission (const vec3d& pos, const vec3d& dir, T_float r, 
				 T_float theta, const T_lambda  em ):
    position (pos), d_ (dir, theta, em), radius_ (r), 
    theta_(theta), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };
  
    virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 
    virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new divergent_floodlight_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** Emission class representing an isotropic external illumination
    field, originating on a sphere and directed inwards.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::external_illumination :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  T_float radius_; ///< The radius of the sphere.
  /** The opening angle of the radiation around the inward direction
      (calculated in the constructor). */
  T_float theta_;
  T_lambda emission_; // The "thing" being emitted.

public:
  /** Creates an external illumination field around the origin. */
  external_illumination (/// Radius where emission originates
			 T_float r, 
			 /** The radius around the origin that should
			     be covered by the emission. */
			 T_float x, 
			 /// The thing being emitted.
			 const T_lambda  em ):
    radius_ (r), theta_(asin(x/r)), emission_ (em) {
    assert(r>=x);
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

    virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, T_float& prob) const; 
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new external_illumination<chromatic_policy, rng_policy_type>(*this));};
};


/** Base class representing volumetric emitters. This class is derived
    from by the concrete classes that implement cells in the various
    grid types. Note that there is a specialization of this class for
    polychromatic runs which can do "effective scattering" but this
    doesn't make sense in the general case since there is no concept
    of "wavelength" unless we're using polychromatic. */
template <typename chromatic_policy, 
	  typename rng_policy_type>
class mcrx::volumetric_emission :
  public isotropic_emission<chromatic_policy, rng_policy_type> {
  friend void assign<>(volumetric_emission&, const volumetric_emission&);
public:
  typedef isotropic_emission<chromatic_policy, rng_policy_type> T_base;
  typedef typename T_base::T_lambda T_lambda;
  typedef typename T_base::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

protected:
  /// The local luminosity (emissivity is this divided by the total mass).
  T_lambda emission_; 

public:
  /** Default constructor creates an empty/invalid emitter. */
  volumetric_emission() : 
    emission_() {};

  volumetric_emission(const T_lambda& em) :
    isotropic_emission<chromatic_policy, rng_policy_type>(em),
    emission_(em) {};

  /* Copy constructor makes an independent copy of the emission array. */
  volumetric_emission(const volumetric_emission& rhs) :
    isotropic_emission<chromatic_policy, rng_policy_type>(rhs),
    emission_(independent_copy(rhs.emission_)) {};

  /** Assignment operator makes a shallow copy. This way we can use
      the assignment operator where we used to use set_data() on the
      grid cell. */
  volumetric_emission& operator=(const volumetric_emission& rhs) {
    isotropic_emission<chromatic_policy, rng_policy_type>::operator=(rhs);
    reference_copy(emission_, rhs.emission_);
    return *this;
  };

  const T_lambda& get_emission() const { return emission_; };
  T_lambda& get_emission() { return emission_; };
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 

  void effective_scatter(T_ray& r, T_lambda& norm, T_rng_policy& rng) const {
    assert(0); };

  /// \name grid_factory stuff. 
  /// These are needed for the grid_factory to work.
  ///@{
  volumetric_emission& operator+= (const volumetric_emission& rhs) {
    emission_+= rhs.emission_; return *this;};
  volumetric_emission& operator*= (const volumetric_emission& rhs) {
    emission_*= rhs.emission_;//content is irrelevant for this operation
    return *this;};
  volumetric_emission& operator/= (T_float rhs) {
    emission_/= rhs;
    return *this;};
  volumetric_emission operator*(const volumetric_emission& rhs) const {
    return volumetric_emission (*this)*= rhs;};
  volumetric_emission operator/(T_float rhs) const {
    return volumetric_emission (*this)/= rhs;};
  ///@}
};


/** Specialization of the volumetric_emission class for polychromatic
    situations. The represent volumetric emitters, typically thermal
    emission from materials. These are distinguished by the fact that
    they are cells in grids, that their emission always is isotropic,
    and that they can perform "effective scattering" where a ray is
    transformed to the shape of the local emissivity while conserving
    energy. This class is derived from by the concrete classes that
    implement cells in the various grid types. */
template <typename rng_policy_type>
class mcrx::volumetric_emission<mcrx::polychromatic_policy, rng_policy_type> :
  public isotropic_emission<polychromatic_policy, rng_policy_type> {
  friend void assign<>(volumetric_emission&, const volumetric_emission&);
public:
  typedef isotropic_emission<polychromatic_policy, rng_policy_type> T_base;
  typedef typename T_base::T_lambda T_lambda;
  typedef typename T_base::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

protected:
  /** Normalizing factor for the emission, such that the integral of
      lnorm_*emission_ = 1. This is only used to avoid having to
      integrate the emission_ over and over for the effective
      scattering. */
  mutable T_float lnorm_;

  /** Pointer to a wavelength vector (this is shared by all emission
      objects in a grid). */
  const T_lambda* lambda_;

  /// The local luminosity (emissivity is this divided by the total mass).
  T_lambda emission_; 

public:
  /** Default constructor creates an empty/invalid emitter. */
  volumetric_emission() : 
    lnorm_(blitz::quiet_NaN(T_float())), lambda_(0), emission_() {};

  volumetric_emission(const T_lambda& em, const T_lambda* lptr) :
    isotropic_emission<polychromatic_policy, rng_policy_type>(em),
    lambda_(lptr), emission_(em) {
    assert(lambda_); emission_update(); };

  /* Copy constructor makes an independent copy of the emission array. */
  volumetric_emission(const volumetric_emission& rhs) :
    isotropic_emission<polychromatic_policy, rng_policy_type>(rhs),
    lnorm_(rhs.lnorm_), lambda_(rhs.lambda_), 
    emission_(independent_copy(rhs.emission_)) {};

  /** Assignment operator makes a shallow copy. This way we can use
      the assignment operator where we used to use set_data() on the
      grid cell. */
  volumetric_emission& operator=(const volumetric_emission& rhs) {
    isotropic_emission<polychromatic_policy, rng_policy_type>::operator=(rhs);
    lnorm_ = rhs.lnorm_;
    lambda_ = rhs.lambda_;
    reference_copy(emission_, rhs.emission_);
    return *this;
  };

  const T_lambda& get_emission() const { return emission_; };
  T_lambda& get_emission() { return emission_; };
  virtual T_float luminosity() const { return 1; };
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 

  /** Signals that the emission (which can reference a master array
      and be updated remotely by e.g. the ir_grid class) has been
      updated and we need to update depedent quantities. */
  void emission_update() { 
    assert(same_size(emission_, *lambda_));
    lnorm_ = 1.0/integrate_quantity(emission_, *lambda_, false);
  };

  /** This function implements "effective scattering", where the
      spectrum of a ray is changed into the source spectrum in the
      cell while conserving energy. The direction is drawn from an
      isotropic distribution. */
  void effective_scatter(T_ray& r, T_lambda& norm, T_rng_policy& rng) const;

  /// \name grid_factory stuff. 
  /// These are needed for the grid_factory to work.
  ///@{
  volumetric_emission& operator+= (const volumetric_emission& rhs) {
    emission_+= rhs.emission_; return *this;};
  volumetric_emission& operator*= (const volumetric_emission& rhs) {
    emission_*= rhs.emission_;//content is irrelevant for this operation
    return *this;};
  volumetric_emission& operator/= (T_float rhs) {
    emission_/= rhs;
    return *this;};
  volumetric_emission operator*(const volumetric_emission& rhs) const {
    return volumetric_emission (*this)*= rhs;};
  volumetric_emission operator/(T_float rhs) const {
    return volumetric_emission (*this)/= rhs;};
  ///@}
};


/** Emission class representing volumetric emission from a rectangular,
    axis-aligned grid cell. IMPORTANT: because this object needs a
    pointer to its cell, we can not simply unify cells and put the
    result in the supercell, as the factory does. After that point,
    the cell_ pointer MUST be restored to point to the current
    cell. Neither the factory nor this object can correct this. It's a
    flaw in the inherent assumption that cell data can always be
    independent of its cell, but this is not true when the position is
    needed.  */
template <typename chromatic_policy, 
	  typename rng_policy_type>
class mcrx::grid_cell_emission :
  public volumetric_emission<chromatic_policy, rng_policy_type> {
  friend void assign<>(grid_cell_emission&, const grid_cell_emission&);
public:
  typedef volumetric_emission<chromatic_policy, rng_policy_type> T_base;
  typedef typename T_base::T_lambda T_lambda;
  typedef typename T_base::T_biaser T_biaser;
  typedef typename T_base::T_ray T_ray;
  typedef rng_policy_type T_rng_policy;

  /** To avoid circularly dependent template parameters, it is
  hardcoded that the grid cells contain this class.  This precludes
  making grids which contain objects which in turn contains
  grid_cell_emission objects, but was chosen for simplicity and the
  lack of the current need for anything more advanced. */
  typedef grid_cell<cell_data<grid_cell_emission, absorber<T_lambda> > > T_cell;

private:
  vec3d min_;
  vec3d size_;

public:
  /** Default constructor creates an empty/invalid emitter. */
  grid_cell_emission() : 
    volumetric_emission<chromatic_policy, rng_policy_type>(),
    min_(0.0), size_(0.0) {};

  /** Constructor makes an emitter with the specified emission in the
      specified cell. Note that emission_ is a reference copy, because
      the ir_grid relies on the individual emitters referencing the
      main L_lambda array. If this is unacceptable, you need to
      explicitly make an independent_copy before passing it to the
      constructor. */
  template<typename T_tracker>
  grid_cell_emission (const T_lambda& em, const T_lambda* lptr,
		      const T_tracker& c):
    volumetric_emission<chromatic_policy, rng_policy_type>(em, lptr),
    min_(c.getmin()), size_(c.getsize()) {};

  /** Constructor makes an emitter with the specified emission in the
      specified cell. This is intended for use with non-polychromatic
      situations where we don't have a concept of wavelength. Note
      that emission_ is a reference copy, because the ir_grid relies
      on the individual emitters referencing the main L_lambda
      array. If this is unacceptable, you need to explicitly make an
      independent_copy before passing it to the constructor. */
  template<typename T_tracker>
  grid_cell_emission (const T_lambda& em, const T_tracker& c):
    volumetric_emission<chromatic_policy, rng_policy_type>(em),
    min_(c.getmin()), size_(c.getsize()) {};

  // copy constructor and assignment operator are implicitly generated

  /** This accessor is so that a grid can use emitters both directly
      and with the the cell_data class, while still retaining that
      possibility. Hopefully use of absorber gets optimized out...  */
  grid_cell_emission& get_emitter() {return *this;};
  const grid_cell_emission& get_emitter() const {return *this;};

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p, 
						T_float& prob) const; 

  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new grid_cell_emission<chromatic_policy, rng_policy_type>(*this));};

  /** Returns the cell pointer. This is needed because if cells are
      unified, we need to be able to restore the cell pointers (very
      hoaky) */
  const T_cell* cell() const { return cell_; };
  /** Updates the emitter with geometry from the specified cell. This
      is needed because if cells are unified, we need to be able to
      restore the cell pointers (very hoaky) */
  template<typename T_tracker>
  void set_cell(const T_tracker& c) { 
    min_ = c.getmin(); size_ = c.getsize(); }

  /** Post_setup is a no-op for this class. (It is necessary because
      of the arepo_cell_emission class.) */
  template<typename T_data>
  static void post_setup(const adaptive_grid<T_data>& g) {};

  /// \name grid_factory stuff. 
  /// These are needed for the grid_factory to work.
  ///@{
  grid_cell_emission& operator+= (const grid_cell_emission& rhs) {
    volumetric_emission<chromatic_policy, rng_policy_type>::operator+=(rhs);
    return *this;};
  grid_cell_emission& operator*= (const grid_cell_emission& rhs) {
    volumetric_emission<chromatic_policy, rng_policy_type>::operator*=(rhs);
    return *this;};
  grid_cell_emission& operator/= (T_float rhs) {
    volumetric_emission<chromatic_policy, rng_policy_type>::operator/=(rhs);
    return *this;};
  grid_cell_emission operator*(const grid_cell_emission& rhs) const {
    return grid_cell_emission (*this)*= rhs;};
  grid_cell_emission operator/(T_float rhs) const {
    return grid_cell_emission (*this)/= rhs;};
  static grid_cell_emission unification (const grid_cell_emission& sum, int n) {
    // emission is extensive
    return sum/n;};
  ///@}
};

template <typename A, typename B>
void mcrx::assign(volumetric_emission<A,B>& lhs,
		  const volumetric_emission<A,B>& rhs) 
{
  assign(lhs.emission_, rhs.emission_);
}  

template <typename A>
void mcrx::assign(volumetric_emission<polychromatic_policy, A>& lhs,
		  const volumetric_emission<polychromatic_policy, A>& rhs) 
{
  lhs.lnorm_ = rhs.lnorm_;
  lhs.lambda_ = rhs.lambda_;
  assign(lhs.emission_, rhs.emission_);
}  

template <typename A, typename B>
void mcrx::assign(grid_cell_emission<A,B>& lhs,
		  const grid_cell_emission<A, B>& rhs) 
{
  assign(static_cast<volumetric_emission<A,B>&>(lhs),
	 static_cast<const volumetric_emission<A,B>&>(rhs));
  lhs.min_=rhs.min_;
  lhs.size_=rhs.size_;
}
  

// ***
// *** function definitions ***
// ***


/** Emits a ray from the particle.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::particle_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::particle_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  ASSERT_ALL(p.data()==p.data());
  ASSERT_ALL(p.data()<1e300);

  // get five random numbers
  blitz::TinyVector<T_float, 5> r;
  rng.rnd(r);

  // draw a radius
  const T_float radius = p.draw_radius(r[0], prob);

  // draw the position in spherical coordinates
  const T_float phi = r [1]*2*constants::pi; 
  const T_float cos_theta = 2*r [2] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);
  prob *= 1/(4*constants::pi);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [3]*2*constants::pi; 
  const T_float dcos_theta = 2*r [4] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);
  const vec3d raydir (dsin_theta*cos (dphi), 
		      dsin_theta*sin (dphi), 
		      dcos_theta);
  // note that there are 2 factors of 1/4pi, one from the position and
  // one from the direction.
  prob *= 1/(4*constants::pi);

  //and create the ray
  boost::shared_ptr<T_ray> rp
    (new T_ray (p.position() + radius*
		vec3d (sin_theta*cos (phi), sin_theta*sin (phi), cos_theta),
		raydir,
		initialize(p.data (), 1.))); // set ray intensity to 0.
  // we do *not* set the doppler factor here because it needs to be
  // emerged to the cameras, etc. instead we return the velocity
  // vector.
  return ray_distribution_pair<T_lambda> (rp, p.data(), this->ad_, p.velocity());
}


/** Isotropically emits a ray from the point source in a random
    direction. The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::pointsource_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::pointsource_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  ASSERT_ALL(emission_==emission_);
  ASSERT_ALL(emission_<1e300);

  // get two random numbers
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  //draw the emission direction (from an isotropic distribution)
  const T_float phi = r [0]*2*constants::pi;  
  const T_float cos_theta = 2*r [1] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);
  prob = 1/(4*constants::pi);

  // and create the ray
  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (position,
		   vec3d (sin_theta*cos (phi),
			  sin_theta*sin (phi),
			  cos_theta),
		   initialize(emission_, 1.))),
      emission_,
      this->ad_);
}


template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::biased_pointsource_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::biased_pointsource_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  // get 2 random numbers
  blitz::TinyVector<T_float, 2> r;
  rng(r);

  //draw the emission direction (NOT from an isotropic distribution)
  const T_float phi = r [0]*2*constants::pi;  
  // this is a P(cos theta) = pi/4*cos([cos theta]*pi/2) distro
  const T_float cos_theta = 2/constants::pi*asin(2*r [1] - 1);
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);
  //const T_float bias = 0.5*pi*cos(cos_theta*0.5*pi);
  const T_float bias = 2/(constants::pi*cos(cos_theta*0.5*constants::pi));

  //prob *= 1/(2*constants::pi)*(constants::pi/4*cos(cos_theta*constants::pi/2));
  // probability needs to be checked if we ever are going to use it
  prob = blitz::quiet_NaN(T_float());

  // create the ray
  boost::shared_ptr<T_ray> ray (new  T_ray (position,
					    vec3d (sin_theta*cos (phi),
						   sin_theta*sin (phi),
						   cos_theta),
					    initialize(emission_, 1.)));
  // the normalization is bias*emission
  return ray_distribution_pair<T_lambda> (ray, bias*emission_, 
					  this->isotropic_ad);
}


/** Emits a ray from the specified point in the specified
    direction. This is trivial because position and direction are
    already specified, all that is needed is to create the ray.  The
    biaser argument is not used.  */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::laser_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::laser_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  prob=1;
  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (position, d.direction (), initialize(emission_, 1.))), 
       emission_, d);
}


/** Emits a ray from in a beam of specified radius in the specified
    direction.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::floodlight_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::floodlight_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  // get 2 random numbers
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  // We need to draw a position on the emitting disk. The radius is
  // drawn from a square-root distribution.
  const T_float rr = sqrt(r [0])*radius;  
  const T_float phi = r [1]*2*constants::pi;
  // Make a position in the xy plane
  const vec3d xypos(rr*cos(phi), rr*sin(phi), 0.);
  // and rotate+translate so the disk is normal to the desired
  // direction and in the desired position
  const vec3d pos = rotate(xypos, d.direction()) + position;

  // probability is 1/area.
  prob = 1.0/(constants::pi*radius*radius);

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos, d.direction (), initialize(emission_, 1.) )), 
       emission_, d);
}


/** Emits a ray from in a beam of specified radius in the specified
    direction.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::divergent_floodlight_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::divergent_floodlight_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  // get 4 random numbers
  blitz::TinyVector<T_float, 4> r;
  rng.rnd(r);

  // We need to draw a position on the emitting disk. The radius is
  // drawn from a square-root distribution.
  const T_float rr = sqrt(r [0])*radius_;  
  const T_float phi = r [1]*2*constants::pi;
  // Make a position in the xy plane
  const vec3d xypos(rr*cos(phi), rr*sin(phi), 0.);
  // and rotate+translate so the disk is normal to the desired
  // direction and in the desired position
  const vec3d pos = rotate(xypos, d_.direction()) + position;

  //draw the emission direction from an isotropic distribution of cos
  //thetas up to theta_, centered around 1 (ie, pointing forward), and
  //rotate into the system of the emission point

  const T_float dphi = r [2]*2*constants::pi;  
  const T_float dcos_theta = 1-(1-cos(theta_))*r [3];
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // probability is 1/area from position and 1/omega from direction.
  const T_float sinhalftheta = sin(0.5*theta_);
  const T_float omega = 4*constants::pi*sinhalftheta*sinhalftheta;
  prob = 1.0/(constants::pi*radius*radius*omega);

  // Make a vector of the emission direction and rotate it into the
  // correct direction
  vec3d dir(dsin_theta*cos (dphi),
	    dsin_theta*sin (dphi),
	    dcos_theta);
  dir=rotate(dir, d_.direction());

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos, dir, initialize(emission_, 1.) )), 
       emission_, d_);
}


/** Emits a ray in an isotropic direction pointing inward from the
    surface of a sphere.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::external_illumination<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::external_illumination<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  // get 4 random numbers
  blitz::TinyVector<T_float, 4> r;
  rng.rnd(r);

  //draw the emission point from an isotropic direction
  const T_float phi = r [0]*2*constants::pi;  
  const T_float cos_theta = 2*r [1] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);

  // Make a vector pointing to the emission point
  const vec3d pos(sin_theta*cos (phi),
		  sin_theta*sin (phi),
		  cos_theta);

  //draw the emission direction from an isotropic distribution of cos
  //thetas up to theta_, centered around -1 (ie, pointing inward), and
  //rotate into the system of the emission point

  const T_float dphi = r [2]*2*constants::pi;  
  const T_float dcos_theta = (1-cos(theta_))*r [3] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // Make a vector of the emission direction and rotate it into the
  // correct direction
  vec3d dir(dsin_theta*cos (dphi),
	    dsin_theta*sin (dphi),
	    dcos_theta);
  dir=rotate(dir, pos);

  // calculate probability if we ever need it
  prob = blitz::quiet_NaN(T_float());

  // generate a cone angular distribution
  const cone_angular_distribution<T_lambda> d(dir, theta_, emission_);

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos*radius_, dir, initialize(emission_, 1.) )), 
       emission_, d);
}


template <typename rng_policy_type>
void
mcrx::volumetric_emission<mcrx::polychromatic_policy, rng_policy_type>::
effective_scatter(T_ray& r, T_lambda& norm, T_rng_policy& rng) const
{
  assert(lnorm_==lnorm_);
  assert(emission_.size()>0);

  // draw direction from an isotropic distribution
  const T_float dphi = rng.rnd()*2*constants::pi; 
  const T_float dcos_theta = 2*rng.rnd() - 1;

  r.effective_scatter(norm, *lambda_, emission_*lnorm_, dcos_theta, dphi);
}


/** Emits a ray from the grid cell. The biaser argument is not used. */
template <typename chromatic_policy,
	  typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::grid_cell_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::grid_cell_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  ASSERT_ALL(this->emission_==this->emission_);
  ASSERT_ALL(this->emission_<1e300);

  // get five random numbers
  blitz::TinyVector<T_float, 5> r;
  rng.rnd(r);
  
  // randomly draw a position
  const vec3d p (r [0], r [1], r [2]);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [3]*2*constants::pi; 
  const T_float dcos_theta = 2*r [4] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // probability is 1/volume from position and 1/4pi from direction.
  prob = 1.0/(product(size_)*4*constants::pi);

  //and create the ray
  return ray_distribution_pair<T_lambda> 
    (boost::shared_ptr<T_ray> 
     (new T_ray (min_ + p*size_,
		 vec3d (dsin_theta*cos (dphi), dsin_theta*sin (dphi), 
			dcos_theta),
		 initialize(this->emission_, 1.))), 
     this->emission_,
     this->ad_);
}

#endif
